/*
 * KineticsSolver.cpp
 *
 *  Created on: 28 Jul 2011
 *      Author: Allan
 */

#include "KineticSolver.h"
#include "Utils/Algorithms.h"
#include "Utils/Assert.h"

KineticSolver::KineticSolver()
{}

KineticSolver::KineticSolver(const Multiphase& multiphase, 
                             const ReactionSystem& eReactions, 
                             const ReactionSystem& kReactions,
                             const vector<string>& jSpecies, 
                             const vector<string>& eSpecies, 
                             const vector<string>& kSpecies)
{
	Initialiaze(multiphase, eReactions, kReactions, jSpecies, eSpecies, kSpecies);
}

KineticSolver::~KineticSolver()
{}

void KineticSolver::Initialiaze(const Multiphase& multiphase, 
                                const ReactionSystem& eReactions, 
                                const ReactionSystem& kReactions,
                                const vector<string>& jSpecies, 
                                const vector<string>& eSpecies, 
                                const vector<string>& kSpecies)
{
	// Set data members
	this->multiphase = multiphase;
	this->eReactions = eReactions;
	this->kReactions = kReactions;
	
	// Set the number of primary (Nj), equilibrium (Ne) and kinetic (Nk) species
	Nj = jSpecies.size();
	Ne = eSpecies.size();
	Nk = kSpecies.size();
	
	// Set the indexes of the primary (jIndexes), equilibrium (eIndexes) and kinetic (kIndexes) species
	jIndexes = multiphase[jSpecies];
	eIndexes = multiphase[eSpecies];
	kIndexes = multiphase[kSpecies];
	
	// Assemble the stoichiometric matrices (uej) and (uee)
	MatrixXd uej = eReactions.AssembleStoichiometricMatrix(jSpecies);
	MatrixXd uee = eReactions.AssembleStoichiometricMatrix(eSpecies);
	
	// Compute (uee_inv), which is the inverse of the stoichiometric matrix (uee)
	FullPivLU<MatrixXd> lu(uee);
	
	MatrixXd uee_inv = lu.inverse();
	
	// Compute the canonical stoichiometric matrix of the equilibrium-controlled reactions
	vej = - uee_inv * uej;
	
	// Initialise the auxiliar Eigen vectors with zero entries
	nj   = VectorXd::Zero(Nj);
	ne   = VectorXd::Zero(Ne);
	nk   = VectorXd::Zero(Nk);
	uj   = VectorXd::Zero(Nj); 
	ujnk = VectorXd::Zero(Nj + Nk);
	
	// Initialiaze the kinetics problem instance 
	problem.Initialiaze(multiphase, eReactions, kReactions, jSpecies, eSpecies, kSpecies);
	
	// Initialise the ODE solver
	odeSolver.Initialiaze(Nj + Nk);
	
	// Initialiaze the equilibrium solver
	eSolver.Initialiaze(multiphase, eReactions, jSpecies, eSpecies);
}

void KineticSolver::SetAlgorithm(OdeAlgorithm algorithm)
{
	odeSolver.Initialiaze(Nj + Nk, algorithm);
}

void KineticSolver::SetControlSolution(double epsAbs, double epsRel)
{
	odeSolver.SetControlSolution(epsAbs, epsRel);
}

bool KineticSolver::Iterate(double& t, double tEnd, double& dt, double T, double P, VectorXd& n)
{
	// Transfer the nj, ne, and nk part of the vector (n) to the vectors (nj), (ne), and (nk) respectively
	for(unsigned j = 0; j < Nj; ++j) nj[j] = n[jIndexes[j]];
	for(unsigned e = 0; e < Ne; ++e) ne[e] = n[eIndexes[e]];
	for(unsigned k = 0; k < Nk; ++k) nk[k] = n[kIndexes[k]];
	
	// Calculate the uj and nk part of the vector (ujnk) = [uj nk], where uj = nj + \sum_e vej * ne
	ujnk.segment(00, Nj) = nj + vej.transpose() * ne;
	ujnk.segment(Nj, Nk) = nk;
	
	// Set the necessary parameters of the kinetic problem object
	problem.SetParameters(T, P, n);
	
	// Use the ODE solver to calculate the vector (ujnk) at the new time (tn) = (t) + (dt)
	odeSolver.Iterate(problem, t, tEnd, dt, ujnk);
	
	// Transfer the uj and nk part of the vector (ujnk) = [uj nk] to the vectors (uj) and (nk) respectively
	uj = ujnk.segment(00, Nj);
	nk = ujnk.segment(Nj, Nk);
	
	// Update the nk part of the vector (n)
	for(unsigned k = 0; k < Nk; ++k) n[kIndexes[k]] = nk[k];
	
	// Solve the equilibrium problem with total composition of the primary species (uj)
	auto diagnostic = eSolver.Solve(uj, T, P, n);
	
	Assert(diagnostic.converged, "Method did not converged.");
	
	return diagnostic.converged;
}

